/*************************************************************************
Purpose: This program calculates the Greenwich phase shifts.
**************************************************************************    
Input: None
**************************************************************************    
Output: None
**************************************************************************    
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

void make_doodson_plots();

void make_doodson_plots(){
  /**********************************************************************
  Purpose:
            This function 
  ***********************************************************************
  Input:  (variable types ending in "_s" are defined in definitions.hpp)
  ***********************************************************************
  Output: (variable types ending in "_s" are defined in definitions.hpp)
  **********************************************************************/
  long long i,j;
  long long interval  = 3600;//In seconds.
  long long numpoints = 48;
  long long offset    = 0;//In seconds.
  
  string x_units_label;
  if(interval == 60) x_units_label = "Time (minutes)";
  else if(interval == 3600) x_units_label = "Time (hours)";
  else if(interval == 3600*24) x_units_label = "Time (days)";
  else if(interval == 3600*24*7) x_units_label = "Time (weeks)";
  else cout<<"WARNING!!! Didn't recognize interval!"<<endl;

  //Needed by write_output().
  results_s temp_results;
  string base = "doodson";
  grid_s temp_grid;
  plot_options_s temp_plot_options;
  ancillary_s temp_ancillary;
  doodson_s doodson;
  
  model_s model;
  model.options.type = 1;
  define_model_titles_and_periods(model);
  
  //Output directly to the output temp folder.
  temp_plot_options.outputfolder = gaiafolder();
  temp_plot_options.outputfolder.append("output/temp/");

  //Override the usual map output.
  temp_results.options.output_choice = 2;

  vector<double> times;
  times.clear();
  //Only a single grid point at a time.
  times.resize(numpoints);
  for(i=0;i<numpoints;i++){
    times.at(i) = (double)(offset + i*interval)/(double)interval;
  }

  //Used to push_back vectors of dbl vectors
  vector<double> temp_v_double;
  //push_back vs of vs of dbl vectors!
  vector< vector<double> > temp_vv_double;
  //push_back vs of string vectors!
  vector<string> temp_v_string;
  
  int plot_counter = -1;//So first plot is at "0".

  /////////////////////////////////////////////////
  //Shailen's Doodson arguments plot.
  /////////////////////////////////////////////////
  plot_counter++;
  doodson.type=1;//Shailen's thesis.
  
  doodson.input_period = m2_period;//Doesn't matter!
  doodson.time = offset;
  greenwich_phase_at_epoch(doodson);//Need this for titles.
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Shailen's Doodson Arguments");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has 6 arguments.
  temp_results.xy.x_values.at(plot_counter).resize(6);
  temp_results.xy.y_values.at(plot_counter).resize(6);

  for(j=1;j<7;j++){
    temp_results.xy.legends.at(plot_counter).push_back(doodson.short_titles.at(j).c_str());
    for(i=0;i<numpoints;i++){
      doodson.time = offset + i*interval;
      greenwich_phase_at_epoch(doodson);
      //They're offset by 1.
      temp_results.xy.x_values.at(plot_counter).at(j-1).push_back(times.at(i)-times.at(0));
      temp_results.xy.y_values.at(plot_counter).at(j-1).push_back(doodson.args.at(j));
    }
  }

  /////////////////////////////////////////////////
  //FES Doodson arguments plot.
  /////////////////////////////////////////////////
  plot_counter++;
  doodson.type=2;//FES code.
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("FES Doodson Arguments");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has 6 arguments.
  temp_results.xy.x_values.at(plot_counter).resize(6);
  temp_results.xy.y_values.at(plot_counter).resize(6);

  for(j=1;j<7;j++){
    temp_results.xy.legends.at(plot_counter).push_back(doodson.short_titles.at(j).c_str());
    for(i=0;i<numpoints;i++){
      doodson.time = offset + i*interval;
      greenwich_phase_at_epoch(doodson);
      //They're offset by 1.
      temp_results.xy.x_values.at(plot_counter).at(j-1).push_back(times.at(i)-times.at(0));
      temp_results.xy.y_values.at(plot_counter).at(j-1).push_back(doodson.args.at(j));
    }
  }
 
  /////////////////////////////////////////////////
  //Simon's Doodson arguments plot.
  /////////////////////////////////////////////////
  plot_counter++;
  doodson.type=3;//Simon 1994.
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Simon's Doodson Arguments");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has 6 arguments.
  temp_results.xy.x_values.at(plot_counter).resize(6);
  temp_results.xy.y_values.at(plot_counter).resize(6);

  for(j=1;j<7;j++){
    temp_results.xy.legends.at(plot_counter).push_back(doodson.short_titles.at(j).c_str());
    for(i=0;i<numpoints;i++){
      doodson.time = offset + i*interval;
      greenwich_phase_at_epoch(doodson);
      //They're offset by 1.
      temp_results.xy.x_values.at(plot_counter).at(j-1).push_back(times.at(i)-times.at(0));
      temp_results.xy.y_values.at(plot_counter).at(j-1).push_back(doodson.args.at(j));
    }
  }

  /////////////////////////////////////////////////
  //Simon - Shailen Doodson arguments plot.
  /////////////////////////////////////////////////
  plot_counter++;
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Simon - Shailen Doodson Arguments");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has 6 arguments.
  temp_results.xy.x_values.at(plot_counter).resize(6);
  temp_results.xy.y_values.at(plot_counter).resize(6);

  for(j=1;j<7;j++){
    temp_results.xy.legends.at(plot_counter).push_back(doodson.short_titles.at(j).c_str());
    for(i=0;i<numpoints;i++){
      //They're offset by 1.
      temp_results.xy.x_values.at(plot_counter).at(j-1).push_back(times.at(i)-times.at(0));
      //-1 is Simon, -3 is Shailen.
      temp_results.xy.y_values.at(plot_counter).at(j-1).push_back(temp_results.xy.y_values.at(plot_counter-1).at(j-1).at(i) - temp_results.xy.y_values.at(plot_counter-3).at(j-1).at(i));
    }
  }

  /////////////////////////////////////////////////
  //Shailen's Greenwich phases plot.
  /////////////////////////////////////////////////
  plot_counter++;
  doodson.type=1;//Shailen's thesis.
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Shailen's Greenwich Phases");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has many timeseries.
  temp_results.xy.x_values.at(plot_counter).resize(model.periods.size());
  temp_results.xy.y_values.at(plot_counter).resize(model.periods.size());

  //Assign legends.
  for(j=0;j<(long long)model.periods.size();j++){
    temp_results.xy.legends.at(plot_counter).push_back(model.titles.at(j).c_str());
  }
  
  //Loop through times.
  for(i=0;i<numpoints;i++){
    doodson.time = offset + i*interval;
    //Loop through periods.
    for(j=0;j<(long long)model.periods.size();j++){
      doodson.input_period = model.periods.at(j);
      greenwich_phase_at_epoch(doodson);
      temp_results.xy.x_values.at(plot_counter).at(j).push_back(times.at(i)-times.at(0));
      temp_results.xy.y_values.at(plot_counter).at(j).push_back(doodson.phase);
    }
  }

  /////////////////////////////////////////////////
  //Simon's Greenwich phases plot.
  /////////////////////////////////////////////////
  plot_counter++;
  doodson.type=3;//Simon 1994.
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Simon's Greenwich Phases");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has many timeseries.
  temp_results.xy.x_values.at(plot_counter).resize(model.periods.size());
  temp_results.xy.y_values.at(plot_counter).resize(model.periods.size());

  //Assign legends.
  for(j=0;j<(long long)model.periods.size();j++){
    temp_results.xy.legends.at(plot_counter).push_back(model.titles.at(j).c_str());
  }
  
  //Loop through times.
  for(i=0;i<numpoints;i++){
    doodson.time = offset + i*interval;
    //Loop through periods.
    for(j=0;j<(long long)model.periods.size();j++){
      doodson.input_period = model.periods.at(j);
      greenwich_phase_at_epoch(doodson);
      temp_results.xy.x_values.at(plot_counter).at(j).push_back(times.at(i)-times.at(0));
      temp_results.xy.y_values.at(plot_counter).at(j).push_back(doodson.phase);
    }
  }


  /////////////////////////////////////////////////
  //Simon - Shailen Doodson arguments plot.
  /////////////////////////////////////////////////
  plot_counter++;
  
  //Loop through all other Doodson numbers because they're all in degrees.
  temp_results.xy.titles.push_back("Simon - Shailen Greenwich Phases");
  temp_results.xy.x_units.push_back(x_units_label);
  temp_results.xy.y_units.push_back("Degrees");

  //This is a brand new plot, so push_back using temp vectors/vvs.
  temp_results.xy.x_values.push_back(temp_vv_double);
  temp_results.xy.y_values.push_back(temp_vv_double);
  temp_results.xy.legends.push_back(temp_v_string);

  //This plot has many timeseries.
  temp_results.xy.x_values.at(plot_counter).resize(model.periods.size());
  temp_results.xy.y_values.at(plot_counter).resize(model.periods.size());

  //Assign legends.
  for(j=0;j<(long long)model.periods.size();j++){
    temp_results.xy.legends.at(plot_counter).push_back(model.titles.at(j).c_str());
  }
  
  //Loop through times.
  for(i=0;i<numpoints;i++){
    //Loop through periods.
    for(j=0;j<(long long)model.periods.size();j++){
      temp_results.xy.x_values.at(plot_counter).at(j).push_back(times.at(i)-times.at(0));
      //-1 is Simon, -2 is Shailen.
      double diff = temp_results.xy.y_values.at(plot_counter-1).at(j).at(i) - temp_results.xy.y_values.at(plot_counter-2).at(j).at(i);
      if(diff>358) diff -= 360;
      if(diff<-358) diff += 360;
      temp_results.xy.y_values.at(plot_counter).at(j).push_back(diff);
    }
  }

  write_output(base, temp_results, temp_grid, temp_plot_options, 1);
}

int main(){
/*
phase diff for   M2:   -6.22, time:  -12 minutes, 52 seconds.
phase diff for   K1:   -9.60, time:  -38 minutes, 16 seconds.
phase diff for   S2:   -0.06, time:  - 6 seconds.
phase diff for   S1:  179.87, time:  11 hours, 59 minutes, 27 seconds.
phase diff for   O1:    3.15, time:  13 minutes, 34 seconds.
phase diff for   P1:   10.16, time:  40 minutes, 44 seconds.
phase diff for   N2:  -63.68, time:  - 2 hours, 14 minutes, 21 seconds.
phase diff for   K2:  160.63, time:   5 hours, 20 minutes, 22 seconds.
phase diff for   Mf:  169.74, time:   6 days, 10 hours, 35 minutes.
phase diff for   Q1:  -54.41, time:  - 4 hours,  3 minutes, 39 seconds.
phase diff for   M4:  -12.16, time:  -12 minutes, 34 seconds.
phase diff for  2N2: -121.52, time:  - 4 hours, 21 minutes, 21 seconds.
phase diff for   Mm:   56.93, time:   4 days,  8 hours, 34 minutes.
phase diff for Msqm:  176.24, time:   3 days, 11 hours, 22 minutes.
phase diff for  Mtm: -132.06, time:  - 3 days,  8 hours, 24 minutes.

New w/o int epoch:

phase diff for   M2: -124.17, time:  - 4 hours, 17 minutes,  2 seconds.
phase diff for   K1:  170.13, time:  11 hours, 18 minutes, 38 seconds.
phase diff for   S2:   -0.06, time:  - 6 seconds.
phase diff for   S1:   -0.13, time:  -32 seconds.
phase diff for   O1:   65.48, time:   4 hours, 41 minutes, 46 seconds.
phase diff for   P1: -169.57, time:  -11 hours, 20 minutes,  7 seconds.
phase diff for   N2:   10.70, time:  22 minutes, 33 seconds.
phase diff for   K2:  160.07, time:   5 hours, 19 minutes, 16 seconds.
phase diff for   Mf:  -72.90, time:  - 2 days, 18 hours, 23 minutes.
phase diff for   Q1: -159.76, time:  -11 hours, 55 minutes, 23 seconds.
phase diff for   M4:  111.98, time:   1 hour, 55 minutes, 54 seconds.
phase diff for  2N2:  145.19, time:   5 hours, 12 minutes, 17 seconds.
phase diff for   Mm: -135.40, time:  - 1 week,  3 days,  8 hours.
phase diff for Msqm:   51.63, time:   1 day, 25 minutes, 19 seconds.
phase diff for  Mtm:  152.97, time:   3 days, 21 hours,  8 minutes.

//*/

  //make_doodson_plots();
  
  time_conversion_s time_conversion; time_conversion.initialized=0;
  time_conversion.calendar_time.resize(6,0);
  
  time_conversion.calendar_time.at(0) = 2005;//Year
  time_conversion.calendar_time.at(1) = 1;//month
  time_conversion.calendar_time.at(2) = 1;//day
  time_conversion.calendar_time.at(3) = 12;//hours
  time_conversion.calendar_time.at(4) = 0;//minutes
  time_conversion.calendar_time.at(5) = 0;//seconds
  cal2utc(time_conversion);//Convert calendar time to UTC seconds.
  //time_conversion.utc_seconds = 0;
  cout<<"UTC seconds: "<<time_conversion.utc_seconds<<endl;

  long long i;
  model_s model;
  model.options.type = 1;
  long long time=time_conversion.utc_seconds;

  define_model_titles_and_periods(model);
  
  cout<<"Time offset: "<<seconds_to_string(time)<<endl;

  astronomicAngle astro;

  double astrotime = (time-fes_utc_int_epoch)*sec2jcy;
  astronomics(astrotime, &astro);
  
  doodson_s doodson;
  doodson.time = time;
  
  /*doodson.type=1;//Shailen's thesis.
  for(i=0;i<(long long)model.periods.size();i++){
    doodson.input_period = model.periods.at(i);
    greenwich_phase_at_epoch(doodson);
    //cout<<"Shailen's phase diff due to beta for        "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<restrict_degrees(doodson.phase)<<" degrees."<<endl;
    cout<<"Shailen's   phase diff due to both for        "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<restrict_degrees(doodson.phase - fes_utc_int_epoch*model.omegas.at(i)*rad2deg)<<" degrees."<<endl;
    //cout<<"Period diff for      "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<doodson.input_period - doodson.output_period<<endl;
    //cout<<"Original period for  "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<seconds_to_string(doodson.input_period)<<endl;
    //cout<<"Shailen's period for "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<seconds_to_string(doodson.output_period)<<endl;
  }//*/
  
  cout<<endl;

  doodson.type=3;//Simon 1994.
  for(i=0;i<(long long)model.periods.size();i++){
    doodson.input_period = model.periods.at(i);
    greenwich_phase_at_epoch(doodson);
    //cout<<"Simon's   phase diff due to beta for        "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<restrict_degrees(doodson.phase)<<" degrees."<<endl;
    //cout<<"Simon's   phase diff due to both for        "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<restrict_degrees(doodson.phase - fes_utc_int_epoch*model.omegas.at(i)*rad2deg)<<" degrees."<<endl;
    //cout<<"Period diff for      "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<doodson.input_period - doodson.output_period<<endl;
    //cout<<"Original period for  "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<seconds_to_string(doodson.input_period)<<endl;
    //cout<<"Simon's period for "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<seconds_to_string(doodson.output_period)<<endl;
    //cout<<"Simon's period for "<<model.titles.at(i)<<": "<<setprecision(10)<<setw(15)<<doodson.output_period<<endl;
  }
  //*/
  cout<<"FES can't calculate phases!"<<endl;
  return 0;
}
